packages = c('tmap', 'SpatialAcc', 'sf', 'ggstatsplot', 'tidyverse')
for(p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p, character.only = T)
}
na_validity_check <-function(input){
print("Check for CRS...")
print(st_crs(input))
print("Check for NA...")
print(input[rowSums(is.na(input))!=0,])
test <- st_is_valid(input,reason=TRUE)
print("Check for invalid geometry...")
length(which(test!= "Valid Geometry"))
}
subzone = st_read(dsn = "data/geospatial/subzone.kml")
## Reading layer `URA_MP19_SUBZONE_NO_SEA_PL' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\subzone.kml' using driver `KML'
## Simple feature collection with 332 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XYZ
## bbox: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
cib = st_read(dsn = "data/geospatial", layer = "cib_gardens")
## Reading layer `cib_gardens' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 1024 features and 14 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 5781.933 ymin: 27790.41 xmax: 45160.86 ymax: 48982.52
## projected CRS: SVY21
comm_club = st_read(dsn = "data/geospatial/community_clubs.kml")
## Reading layer `COMMUNITYCLUBS' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\community_clubs.kml' using driver `KML'
## Simple feature collection with 119 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 103.6923 ymin: 1.274863 xmax: 103.9592 ymax: 1.451682
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
comm_use_site = st_read(dsn = "data/geospatial/community_use_sites.kml")
## Reading layer `COMM_USE_SITES' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\community_use_sites.kml' using driver `KML'
## Simple feature collection with 250 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XYZ
## bbox: xmin: 103.6228 ymin: 1.269832 xmax: 103.9653 ymax: 1.461749
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
dus_sports_facilities = st_read(dsn = "data/geospatial", layer = "dus_school_sports_facilities")
## Reading layer `dus_school_sports_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 174 features and 15 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11739.27 ymin: 28604.3 xmax: 42410.51 ymax: 48696.35
## projected CRS: SVY21
nature_areas = st_read(dsn = "data/geospatial/nature_areas.kml")
## Reading layer `UP_G_MP19_PKWB_NATURE_PL' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nature_areas.kml' using driver `KML'
## Simple feature collection with 62 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XYZ
## bbox: xmin: 103.644 ymin: 1.247164 xmax: 104.0767 ymax: 1.451329
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
np_activity_area = st_read(dsn = "data/geospatial/nparks_activity_area.kml")
## Reading layer `NPARKS_ACTIVITY_AREA' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nparks_activity_area.kml' using driver `KML'
## Simple feature collection with 599 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XYZ
## bbox: xmin: 103.693 ymin: 1.262527 xmax: 103.9936 ymax: 1.462811
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
np_parks = st_read(dsn = "data/geospatial/nparks_parks.kml")
## Reading layer `NPARKS_PARKS' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nparks_parks.kml' using driver `KML'
## Simple feature collection with 390 features and 2 fields
## geometry type: MULTIPOLYGON
## dimension: XYZ
## bbox: xmin: 103.6925 ymin: 1.265889 xmax: 104.008 ymax: 1.46419
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
np_play_fitness = st_read(dsn = "data/geospatial/nparks_play_fitness_equipment.kml")
## Reading layer `NPARKS_PLAY_FITNESS_EQ' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\nparks_play_fitness_equipment.kml' using driver `KML'
## Simple feature collection with 3108 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 103.6929 ymin: 1.262622 xmax: 103.9935 ymax: 1.457533
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
park_loop = st_read(dsn = "data/geospatial/park_connector_loop.kml")
## Reading layer `PARK_CONNECTOR_LOOP' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\park_connector_loop.kml' using driver `KML'
## Simple feature collection with 233 features and 2 fields
## geometry type: LINESTRING
## dimension: XYZ
## bbox: xmin: 103.6933 ymin: 1.272949 xmax: 104.0316 ymax: 1.462852
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
preschools = st_read(dsn = "data/geospatial/preschools.kml")
## Reading layer `PRESCHOOLS_LOCATION' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\preschools.kml' using driver `KML'
## Simple feature collection with 1925 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 103.6824 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
schools = st_read(dsn = "data/geospatial", layer = "schools")
## Reading layer `schools' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 327 features and 43 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11799.7 ymin: 28571.39 xmax: 42434.9 ymax: 48727.31
## projected CRS: SVY21 / Singapore TM
student_care = st_read(dsn = "data/geospatial/student_care.kml")
## Reading layer `STUDENTCARE' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\student_care.kml' using driver `KML'
## Simple feature collection with 425 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 103.6877 ymin: 1.27474 xmax: 103.963 ymax: 1.456667
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
sportsg_facilities= st_read(dsn = "data/geospatial/sportsg_sports_facilities.kml")
## Reading layer `SPORTSG_SPORT_FACILITIES' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\sportsg_sports_facilities.kml' using driver `KML'
## Simple feature collection with 35 features and 2 fields
## geometry type: POLYGON
## dimension: XYZ
## bbox: xmin: 103.6932 ymin: 1.285977 xmax: 103.9526 ymax: 1.435855
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
watersports_facilities = st_read(dsn = "data/geospatial", layer = "water_sport_facilities")
## Reading layer `water_sport_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 17 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 12469.6 ymin: 30191.07 xmax: 45423.18 ymax: 48969.16
## projected CRS: SVY21
na_validity_check(subzone)
## [1] "Check for CRS..."
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
## [1] "Check for NA..."
## Simple feature collection with 0 features and 2 fields
## bbox: xmin: NA ymin: NA xmax: NA ymax: NA
## geographic CRS: WGS 84
## [1] Name Description geometry
## <0 rows> (or 0-length row.names)
## [1] "Check for invalid geometry..."
## [1] 6
subzone <- st_set_crs(subzone, 3414)
subzone<- st_make_valid(subzone)
getName <- function(input){
noHtml <- str_replace_all(input,"<.*?>","")
matched_str <- str_match(noHtml, "\\sSUBZONE_N (.*?) SUBZONE_C")[,2]
return(str_trim(matched_str))
}
subzone <- subzone %>%
mutate(`SUBZONE_N` = getName(subzone$Description))
head(subzone, 1)
nature_areas <- st_transform(nature_areas, 3414)
nature_areas<- st_make_valid(nature_areas)
np_activity_area <- st_transform(np_activity_area, 3414)
np_parks <- st_transform(np_parks, 3414)
np_parks<- st_make_valid(np_parks)
np_play_fitness <- st_transform(np_play_fitness, 3414)
park_loop <- st_transform(park_loop, 3414)
preschools <- st_transform(preschools, 3414)
sportsg_facilities <- st_transform(sportsg_facilities, 3414)
student_care <- st_transform(student_care, 3414)
watersports_facilities <- st_set_crs(watersports_facilities, 3414)
cib <- st_set_crs(cib, 3414)
comm_club <- st_transform(comm_club, 3414)
comm_use_site <- st_transform(comm_use_site, 3414)
dus_sports_facilities <- st_set_crs(dus_sports_facilities, 3414)
HDB
hdb <- read_csv('data/aspatial/hdb_osm.csv')
## Parsed with column specification:
## cols(
## .default = col_logical(),
## X = col_double(),
## Y = col_double(),
## full_id = col_character(),
## osm_id = col_double(),
## osm_type = col_character(),
## addr_city = col_character(),
## addr_house = col_character(),
## addr_postc = col_character(),
## addr_stree = col_character(),
## building = col_character(),
## building_l = col_character(),
## name = col_character(),
## type = col_character(),
## addr_count = col_character(),
## residentia = col_character(),
## year = col_double(),
## addr_hou_1 = col_character(),
## building_c = col_character(),
## building_u = col_character(),
## building_m = col_character()
## # ... with 41 more columns
## )
## See spec(...) for full column specifications.
## Warning: 819 parsing failures.
## row col expected actual file
## 1029 amenity 1/0/T/F/TRUE/FALSE parking 'data/aspatial/hdb_osm.csv'
## 1029 parking 1/0/T/F/TRUE/FALSE multi-storey 'data/aspatial/hdb_osm.csv'
## 1193 operator_2 1/0/T/F/TRUE/FALSE HDB 'data/aspatial/hdb_osm.csv'
## 1194 operator_2 1/0/T/F/TRUE/FALSE HDB 'data/aspatial/hdb_osm.csv'
## 1195 operator_2 1/0/T/F/TRUE/FALSE HDB 'data/aspatial/hdb_osm.csv'
## .... .......... .................. ............ ...........................
## See problems(...) for more details.
hdb <- st_as_sf(hdb,
coords = c('X', 'Y'),
crs = 4326) %>%
st_transform(3414)
population <- read_csv('data/aspatial/population.csv')
## Parsed with column specification:
## cols(
## PA = col_character(),
## SZ = col_character(),
## AG = col_character(),
## Sex = col_character(),
## TOD = col_character(),
## Pop = col_double(),
## Time = col_double()
## )
\[Number\ of\ children\ aged\ 0\ to\ 7\ in\ subzone\ =\ Number\ of\ children\ aged\ 0\ to\ 4 + (\frac{3}{5}*Number\ of\ children\ aged\ 5\ to\ 9)\]
children_population <- population %>%
filter(Time == 2020) %>%
group_by(PA, SZ, AG, TOD) %>%
summarise(`POP` = sum(`Pop`)) %>%
ungroup() %>%
pivot_wider(names_from = AG, values_from = POP) %>%
mutate(children_population = `0_to_4` + (3/5 * `5_to_9`)) %>%
select(PA, SZ, TOD, children_population) %>%
pivot_wider(names_from = TOD, values_from = children_population) %>%
mutate_at(.vars = vars(PA, SZ), .funs = ~ toupper(.)) %>%
mutate(children_hdb = rowSums(.[4:7])) %>%
select(SZ, children_hdb)
## `summarise()` regrouping output by 'PA', 'SZ', 'AG' (override with `.groups` argument)
children_population
Join population data with subzone data
subzone_children <- subzone %>%
left_join(children_population, by = c('SUBZONE_N' = 'SZ')) %>%
select(-c(Description, Name))
glimpse(subzone_children)
## Rows: 332
## Columns: 3
## $ SUBZONE_N <chr> "MARINA EAST", "INSTITUTION HILL", "ROBERTSON QUAY", "...
## $ children_hdb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 492, 0, 358,...
## $ geometry <GEOMETRY [m]> MULTIPOLYGON (((103.8802 1...., MULTIPOLYGON ...
Calculate number of children aged 0 to 7 in each HDB block
\[Number\ of\ children\ in\ each\ HDB\ =\ \frac{Number\ of\ children\ in\ subzone}{Number\ of\ HDB\ blocks\ in\ subzone}\]
Count number of HDB blocks in each subzone
subzone_children$num_hdb <- lengths(st_intersects(subzone_children, hdb))
summary(subzone_children$num_hdb)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
qtm(subzone_children, 'num_hdb')
qtm(subzone_children, 'children_hdb')
Calculate number of children in each HDB: demand
subzone_children <- subzone_children %>%
# Calculate demand. If there are no HDB
mutate(children_each_hdb = ifelse(num_hdb != 0, children_hdb / num_hdb, 0))
summary(subzone_children$children_each_hdb)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
tm_shape(subzone_children) +
tm_polygons('children_each_hdb',
style = 'jenks')
Add information on number of children for each HDB to HDB data
hdb <- hdb %>%
select(`ID`, `addr_house`, `addr_postc`, `addr_stree`, `building_l`) %>%
rename(house_number = addr_house,
postal_code = addr_postc,
street = addr_stree,
num_levels = building_l) %>%
st_join(subzone_children, join = st_intersects)
glimpse(hdb)
## Rows: 9,806
## Columns: 10
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15...
## $ house_number <chr> "21", "21", "94B", "94B", "94C", "94C", "25B", "2...
## $ postal_code <chr> "270021", "270021", "461094", "461094", "462094",...
## $ street <chr> "Ghim Moh Road", "Ghim Moh Road", "Bedok North Av...
## $ num_levels <chr> "10", "10", "13", "13", "13", "13", "27", "27", N...
## $ SUBZONE_N <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ children_hdb <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ num_hdb <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ children_each_hdb <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ geometry <POINT [m]> POINT (23015.11 32522.95), POINT (23015.11 ...
student_care <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'student_care') %>%
select(ID, Name, ADDRESSSTR, ADDRESSPOS) %>%
rename(name = Name,
address = ADDRESSSTR,
postal_code = ADDRESSPOS)
## Reading layer `student_care' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 425 features and 28 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11796.89 ymin: 28579.85 xmax: 42433.39 ymax: 48696.35
## projected CRS: SVY21 / Singapore TM
glimpse(student_care)
## Rows: 425
## Columns: 5
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name <chr> "Ananias Centre (Choa Chu Kang)", "Ananias Centre (Clem...
## $ address <chr> "Blk 289 Choa Chu Kang Avenue 3 #01-262 SINGAPORE 68028...
## $ postal_code <chr> "680289", "121440", "760160", "760417", "530540", "6801...
## $ geometry <POINT [m]> POINT (17694.75 40144.54), POINT (20305.83 33142....
dm_student_care <- read_csv('data/aspatial/distance matrix/hdb_studentcare.csv')
## Parsed with column specification:
## cols(
## origin_id = col_double(),
## destination_id = col_double(),
## entry_cost = col_double(),
## network_cost = col_double(),
## exit_cost = col_double(),
## total_cost = col_double()
## )
glimpse(dm_student_care)
## Rows: 4,167,550
## Columns: 6
## $ origin_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost <dbl> 13126.843, 13126.843, 19604.224, 19672.647, 17050.73...
## $ exit_cost <dbl> 35.10268, 35.10268, 36.57691, 98.50361, 146.68192, 3...
## $ total_cost <dbl> 13251.237, 13251.237, 19730.092, 19860.441, 17286.70...
Tidy matrix – make wider
dm_student_care <- dm_student_care %>%
# Obtain final matrix
select(origin_id, destination_id, total_cost) %>%
pivot_wider(names_from = destination_id, values_from = total_cost) %>%
select(-origin_id)
# Convert to km
dm_student_care <- as.matrix(dm_student_care/1000)
# dm_student_care
Calculate average student care capacity for each student care centre
\[Average\ student\ care\ capacity\ for\ each\ centre\ = \frac{Total\ capacity\ of\ student\ care\ in\ Singapore}{Number\ of\ student\ care\ centres\ in\ Singapore}\]
student_care_capacity <- round(42907 / nrow(student_care), 0)
student_care_capacity
## [1] 101
Calculate Hansen accessibility
student_care_capacity_list <- rep(student_care_capacity, nrow(student_care))
acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
student_care_capacity_list,
dm_student_care,
d0 = 50,
power = 1.5,
family = "Hansen"))
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_student_care <- bind_cols(hdb, acc_Hansen)
hansen_student_care
Import schools
pri_schools <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'schools_primary') %>%
select(ID, school_nam, address, postal_cod) %>%
rename(name = school_nam,
postal_code = postal_cod)
## Reading layer `schools_primary' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 170 features and 44 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11799.7 ymin: 28571.39 xmax: 42434.9 ymax: 48727.31
## projected CRS: SVY21 / Singapore TM
glimpse(pri_schools)
## Rows: 170
## Columns: 5
## $ ID <dbl> 1, 3, 5, 6, 7, 10, 16, 17, 18, 22, 24, 28, 30, 31, 37, ...
## $ name <chr> "ADMIRALTY PRIMARY SCHOOL", "AHMAD IBRAHIM PRIMARY SCHO...
## $ address <chr> "11 WOODLANDS CIRCLE", "10 YISHUN STREET 11", "100 ...
## $ postal_code <chr> "738907", "768643", "579646", "159016", "544969", "5699...
## $ geometry <POINT [m]> POINT (24330.56 47178.67), POINT (27932.25 46173....
Import distance matrix
dm_pri_schools <- read_csv('data/aspatial/distance matrix/hdb_school.csv')
## Parsed with column specification:
## cols(
## origin_id = col_double(),
## destination_id = col_double(),
## entry_cost = col_double(),
## network_cost = col_double(),
## exit_cost = col_double(),
## total_cost = col_double()
## )
glimpse(dm_pri_schools)
## Rows: 1,667,020
## Columns: 6
## $ origin_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 3, 5, 6, 7, 10, 16, 17, 18, 22, 24, 28, 30, 31, 3...
## $ entry_cost <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost <dbl> 19430.115, 19880.538, 11378.661, 6226.051, 18680.708...
## $ exit_cost <dbl> 88.249817, 56.515963, 80.705633, 7.244451, 93.627581...
## $ total_cost <dbl> 19607.656, 20026.345, 11548.657, 6322.587, 18863.627...
Tidy matrix – make wider
dm_pri_schools <- dm_pri_schools %>%
# Obtain final matrix
select(origin_id, destination_id, total_cost) %>%
pivot_wider(names_from = destination_id, values_from = total_cost) %>%
select(-origin_id)
# Convert to km
dm_pri_schools <- as.matrix(dm_pri_schools/1000)
Calculate average primary school capacity for each school
school_capacity <- round(37671 / nrow(pri_schools))
school_capacity
## [1] 222
Calculate Hansen accessibility
school_capacity_list <- rep(school_capacity, nrow(pri_schools))
acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
school_capacity_list,
dm_pri_schools,
d0 = 50,
power = 2,
family = "Hansen"))
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_pri_schools <- bind_cols(hdb, acc_Hansen)
hansen_pri_schools
Import water sports facilities
watersports_facilities <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'water_sport_facilities') %>%
select(ID, NAME, ADDRESSSTR, ADDRESSPOS) %>%
rename(name = NAME,
address = ADDRESSSTR,
postal_code = ADDRESSPOS)
## Reading layer `water_sport_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 12469.6 ymin: 30191.07 xmax: 45423.18 ymax: 48969.16
## projected CRS: SVY21 / Singapore TM
glimpse(watersports_facilities)
## Rows: 32
## Columns: 5
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name <chr> "Toa Payoh Sports and Recreation Centre", "Water-Ventur...
## $ address <chr> "Toa Payoh Lor. 6", "East Coast Parkway", "Bishan Stree...
## $ postal_code <int> 319392, 468961, 579778, 648965, 339627, 436752, 148948,...
## $ geometry <POINT [m]> POINT (29869.2 34746.17), POINT (41233.04 32670.0...
Import distance matrix
dm_watersports_facilities <- read_csv('data/aspatial/distance matrix/hdb_watersports.csv')
## Parsed with column specification:
## cols(
## origin_id = col_double(),
## destination_id = col_double(),
## entry_cost = col_double(),
## network_cost = col_double(),
## exit_cost = col_double(),
## total_cost = col_double()
## )
glimpse(dm_watersports_facilities)
## Rows: 313,792
## Columns: 6
## $ origin_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost <dbl> 9060.498, 20286.510, 10517.724, 10517.724, 11776.135...
## $ exit_cost <dbl> 57.776199, 38.678467, 86.974437, 86.974437, 31.99397...
## $ total_cost <dbl> 9207.565, 20414.480, 10693.990, 10693.990, 11897.420...
Tidy matrix – make wider
dm_watersports_facilities <- dm_watersports_facilities %>%
# Obtain final matrix
select(origin_id, destination_id, total_cost) %>%
pivot_wider(names_from = destination_id, values_from = total_cost) %>%
select(-origin_id)
# Convert to km
dm_watersports_facilities <- as.matrix(dm_watersports_facilities/1000)
Calculate Hansen accessibility
watersports_capacity_list <- rep(1, nrow(watersports_facilities))
acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
watersports_capacity_list,
dm_watersports_facilities,
d0 = 50,
power = 1.5,
family = "Hansen"))
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_watersports_facilities <- bind_cols(hdb, acc_Hansen)
hansen_watersports_facilities
Import DUS sports facilities
dus_sports_facilities <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'dus_school_sports_facilities') %>%
select(ID, SCHOOL_NAM, ADDRESS, POSTAL_COD, FACILITIES) %>%
rename(name = SCHOOL_NAM,
address = ADDRESS,
postal_code = POSTAL_COD,
facilities = FACILITIES)
## Reading layer `dus_school_sports_facilities' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 174 features and 16 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11739.27 ymin: 28604.3 xmax: 42410.51 ymax: 48696.35
## projected CRS: SVY21 / Singapore TM
glimpse(dus_sports_facilities)
## Rows: 174
## Columns: 6
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name <chr> "Changkat Changi Secondary School", "Chestnut Drive Sec...
## $ address <chr> "23 Simei Street 3", "58 Chestnut Drive", "20 Choa Chu ...
## $ postal_code <chr> "529894", "679301", "689905", "688845", "129903", "1299...
## $ facilities <chr> "Chargeable Field", "Chargeable Field, Basketball Court...
## $ geometry <POINT [m]> POINT (41288.67 35836.16), POINT (21042.42 38992....
Import distance matrix
dm_dus_school_sports_facilities <- read_csv('data/aspatial/distance matrix/hdb_dus.csv')
## Parsed with column specification:
## cols(
## origin_id = col_double(),
## destination_id = col_double(),
## entry_cost = col_double(),
## network_cost = col_double(),
## exit_cost = col_double(),
## total_cost = col_double()
## )
glimpse(dm_dus_school_sports_facilities)
## Rows: 1,706,244
## Columns: 6
## $ origin_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost <dbl> 22496.583, 9295.449, 13681.212, 11848.690, 11848.690...
## $ exit_cost <dbl> 100.22417, 45.43876, 80.20923, 39.46815, 39.46815, 1...
## $ total_cost <dbl> 22686.099, 9430.179, 13850.712, 11977.449, 11977.449...
Tidy matrix – make wider
dm_dus_school_sports_facilities <- dm_dus_school_sports_facilities %>%
# Obtain final matrix
select(origin_id, destination_id, total_cost) %>%
pivot_wider(names_from = destination_id, values_from = total_cost) %>%
select(-origin_id)
# Convert to km
dm_dus_school_sports_facilities <- as.matrix(dm_dus_school_sports_facilities/1000)
Calculate Hansen accessibility
dus_capacity_list <- rep(1, nrow(dus_sports_facilities))
acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
dus_capacity_list,
dm_dus_school_sports_facilities,
d0 = 50,
power = 2,
family = "Hansen"))
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_dus_school_sports_facilities <- bind_cols(hdb, acc_Hansen)
hansen_dus_school_sports_facilities
Import community in bloom gardens
cib_gardens <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'cib_gardens') %>%
select(ID, ADDRESS, DIVISION, CONSTITUEN, CATEGORY) %>%
rename(address = ADDRESS,
division = DIVISION,
constituency = CONSTITUEN,
category = CATEGORY)
## Reading layer `cib_gardens' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 1024 features and 15 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 5781.933 ymin: 27790.41 xmax: 45160.86 ymax: 48982.52
## projected CRS: SVY21 / Singapore TM
glimpse(cib_gardens)
## Rows: 1,024
## Columns: 6
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,...
## $ address <chr> "Blk 730 Tampines St 71 #01-61", "TAI HWAN TERRACE, BL...
## $ division <chr> "TAMPINES CENTRAL", "SERANGOON", "TAMPINES CENTRAL", "...
## $ constituency <chr> "TAMPINES GRC", "ALJUNIED GRC", "TAMPINES GRC", "TAMPI...
## $ category <chr> "Public Estate", "Private Estate", "Public Estate", "P...
## $ geometry <POINT [m]> POINT (39151.3 37639.3), POINT (31058.07 37624.8...
Import distance matrix
dm_cib_gardens <- read_csv('data/aspatial/distance matrix/hdb_cib.csv')
## Parsed with column specification:
## cols(
## origin_id = col_double(),
## destination_id = col_double(),
## entry_cost = col_double(),
## network_cost = col_double(),
## exit_cost = col_double(),
## total_cost = col_double()
## )
glimpse(dm_cib_gardens)
## Rows: 10,041,344
## Columns: 6
## $ origin_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost <dbl> 89.29121, 89.29121, 89.29121, 89.29121, 89.29121, 89...
## $ network_cost <dbl> 19708.18, 13155.52, 19763.04, 21985.90, 21768.26, 22...
## $ exit_cost <dbl> 44.587488, 4.749120, 109.933316, 16.257711, 67.14973...
## $ total_cost <dbl> 19842.06, 13249.56, 19962.26, 22091.45, 21924.70, 22...
Tidy matrix – make wider
dm_cib_gardens <- dm_cib_gardens %>%
# Obtain final matrix
select(origin_id, destination_id, total_cost) %>%
pivot_wider(names_from = destination_id, values_from = total_cost) %>%
select(-origin_id)
# Convert to km
dm_cib_gardens <- as.matrix(dm_cib_gardens/1000)
Calculate Hansen accessibility
cib_capacity_list <- rep(1, nrow(cib_gardens))
acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
cib_capacity_list,
dm_cib_gardens,
d0 = 50,
power = 2.5,
family = "Hansen"))
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_cib_gardens <- bind_cols(hdb, acc_Hansen)
hansen_cib_gardens
Import pre-schools
preschools <- st_read(dsn = 'data/geospatial/distance_matrix_', layer = 'preschools') %>%
select(ID, CENTRE_NAM, ADDRESS, POSTAL_COD) %>%
rename(name = CENTRE_NAM,
address = ADDRESS,
postal_code = POSTAL_COD)
## Reading layer `preschools' from data source `C:\Users\Xiao Rong\Desktop\School\Geospatial Analytics and Applications\Project\gis-project\analysis\data\geospatial\distance_matrix_' using driver `ESRI Shapefile'
## Simple feature collection with 1925 features and 19 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11203.01 ymin: 25596.33 xmax: 45404.24 ymax: 49300.88
## projected CRS: SVY21 / Singapore TM
glimpse(preschools)
## Rows: 1,925
## Columns: 5
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ name <chr> "BRILLIANT TOTS PTE. LTD.", "BUBBLESLAND PLAYHOUSE PTE ...
## $ address <chr> "610, JURONG WEST STREET 65, #01 - 534, S 640610", "238...
## $ postal_code <chr> "640610", "540238", "737856", "730369", "542327", "8212...
## $ geometry <POINT [m]> POINT (13258.34 35611.04), POINT (35272.09 41373....
Import distance matrix
dm_preschools <- read_csv('data/aspatial/distance matrix/hdb_preschools.csv')
## Parsed with column specification:
## cols(
## origin_id = col_double(),
## destination_id = col_double(),
## entry_cost = col_double(),
## network_cost = col_double(),
## exit_cost = col_double(),
## total_cost = col_double()
## )
glimpse(dm_preschools)
## Rows: 18,876,550
## Columns: 6
## $ origin_id <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ destination_id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ entry_cost <dbl> NA, 89.29121, 89.29121, 89.29121, 89.29121, 89.29121...
## $ network_cost <dbl> NA, 19008.163, 19731.005, 17167.844, 19267.844, 1926...
## $ exit_cost <dbl> NA, 57.89667, 81.98225, 72.21746, 116.03101, 116.031...
## $ total_cost <dbl> NA, 19155.351, 19902.279, 17329.353, 19473.166, 1947...
Tidy matrix – make wider
dm_preschools <- dm_preschools %>%
# Obtain final matrix
select(origin_id, destination_id, total_cost) %>%
pivot_wider(names_from = destination_id, values_from = total_cost) %>%
select(-origin_id)
# Convert to km
dm_preschools <- as.matrix(dm_preschools/1000)
Calculate average pre-school capacity
preschool_capacity <- round(215579 / nrow(preschools))
preschool_capacity
## [1] 112
Calculate Hansen accessibility
preschool_capacity_list <- rep(preschool_capacity, nrow(preschools))
acc_Hansen <- data.frame(ac(hdb$children_each_hdb,
preschool_capacity_list,
dm_preschools,
d0 = 50,
power = 2,
family = "Hansen"))
colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- as_tibble(acc_Hansen)
hansen_preschools <- bind_cols(hdb, acc_Hansen)
hansen_preschools
hansen_preschools[is.na(hansen_preschools)] <- 0
glimpse(hansen_preschools)
## Rows: 9,806
## Columns: 11
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15...
## $ house_number <chr> "21", "21", "94B", "94B", "94C", "94C", "25B", "2...
## $ postal_code <chr> "270021", "270021", "461094", "461094", "462094",...
## $ street <chr> "Ghim Moh Road", "Ghim Moh Road", "Bedok North Av...
## $ num_levels <chr> "10", "10", "13", "13", "13", "13", "27", "27", "...
## $ SUBZONE_N <chr> "0", "0", "0", "0", "0", "0", "0", "0", "0", "0",...
## $ children_hdb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ num_hdb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ children_each_hdb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ accHansen <dbl> 0.00000, 316.59300, 36.35925, 36.56049, 191.26065...
## $ geometry <POINT [m]> POINT (23015.11 32522.95), POINT (23015.11 ...
Function to compute risk weights
compute_risk_weights <- function(df, name) {
first_quantile_hansen <- quantile(df$accHansen, 0.25)[[1]]
median_hansen <- median(df$accHansen)
third_quantile_hansen <- quantile(df$accHansen, 0.75)[[1]]
max_hansen <- max(df$accHansen)
df1 <- df %>%
mutate('risk_{name}' := case_when(
accHansen < first_quantile_hansen ~ 4,
(accHansen >= first_quantile_hansen) & (accHansen < median_hansen) ~ 3,
(accHansen >= median_hansen) & (accHansen < third_quantile_hansen) ~ 2,
(accHansen >= third_quantile_hansen) ~ 1))
return(df1)
}
Compute risk weights
hansen_cib_gardens <- compute_risk_weights(hansen_cib_gardens, 'cib_gardens')
hansen_dus_school_sports_facilities <- compute_risk_weights(hansen_dus_school_sports_facilities, 'dus_sports')
hansen_preschools <- compute_risk_weights(hansen_preschools, 'preschools')
hansen_pri_schools <- compute_risk_weights(hansen_pri_schools, 'pri_schools')
hansen_student_care <- compute_risk_weights(hansen_student_care, 'student_care')
hansen_watersports_facilities <- compute_risk_weights(hansen_watersports_facilities, 'watersports_facilities')
hdb_risk <- hdb %>%
cbind(hansen_cib_gardens$risk_cib_gardens,
hansen_dus_school_sports_facilities$risk_dus_sports,
hansen_preschools$risk_preschools,
hansen_pri_schools$risk_pri_schools,
hansen_student_care$risk_student_care,
hansen_watersports_facilities$risk_watersports_facilities) %>%
rename(risk_cib_gardens = hansen_cib_gardens.risk_cib_gardens,
risk_dus_sports = hansen_dus_school_sports_facilities.risk_dus_sports,
risk_preschools = hansen_preschools.risk_preschools,
risk_pri_schools = hansen_pri_schools.risk_pri_schools,
risk_student_care = hansen_student_care.risk_student_care,
risk_watersports_facilities = hansen_watersports_facilities.risk_watersports_facilities) %>%
mutate(physical_health_risk = (risk_dus_sports + risk_watersports_facilities) / 2) %>%
mutate(social_competence_risk = (risk_dus_sports + risk_preschools + risk_pri_schools + risk_student_care + risk_watersports_facilities) / 5) %>%
mutate(emotional_maturity_risk = risk_cib_gardens) %>%
mutate(composite_risk = (physical_health_risk + social_competence_risk + emotional_maturity_risk) / 3)
glimpse(hdb_risk)
## Rows: 9,806
## Columns: 20
## $ ID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...
## $ house_number <chr> "21", "21", "94B", "94B", "94C", "94C",...
## $ postal_code <chr> "270021", "270021", "461094", "461094",...
## $ street <chr> "Ghim Moh Road", "Ghim Moh Road", "Bedo...
## $ num_levels <chr> "10", "10", "13", "13", "13", "13", "27...
## $ SUBZONE_N <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ children_hdb <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ num_hdb <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ children_each_hdb <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
## $ risk_cib_gardens <dbl> 3, 3, 3, 3, 2, 2, 1, 1, 2, 2, 2, 2, 2, ...
## $ risk_dus_sports <dbl> 3, 3, 3, 3, 2, 2, 3, 3, 3, 3, 1, 1, 1, ...
## $ risk_preschools <dbl> 4, 3, 4, 4, 4, 4, 3, 3, 2, 3, 3, 2, 2, ...
## $ risk_pri_schools <dbl> 2, 2, 3, 3, 3, 3, 4, 4, 1, 2, 2, 1, 1, ...
## $ risk_student_care <dbl> 3, 3, 2, 2, 3, 3, 3, 3, 2, 2, 2, 1, 1, ...
## $ risk_watersports_facilities <dbl> 1, 1, 2, 2, 2, 2, 4, 4, 1, 1, 1, 1, 1, ...
## $ geometry <POINT [m]> POINT (23015.11 32522.95), POINT ...
## $ physical_health_risk <dbl> 2.0, 2.0, 2.5, 2.5, 2.0, 2.0, 3.5, 3.5,...
## $ social_competence_risk <dbl> 2.6, 2.4, 2.8, 2.8, 2.8, 2.8, 3.4, 3.4,...
## $ emotional_maturity_risk <dbl> 3, 3, 3, 3, 2, 2, 1, 1, 2, 2, 2, 2, 2, ...
## $ composite_risk <dbl> 2.533333, 2.466667, 2.766667, 2.766667,...
tm_shape(subzone_children) +
tm_polygons(col = 'gray90') +
tm_shape(hdb_risk) +
tm_bubbles(col = 'composite_risk',
size = 0.1,
alpha = 0.5,
border.lwd = NA)
# st_write(hdb_risk, '../shinyapp/data/geospatial/hdb_risk.shp', driver = 'ESRI Shapefile')
# st_write(subzone_children, '../shinyapp/data/geospatial/subzone_children.shp', driver = 'ESRI Shapefile')
tmap_mode('view')
## tmap mode set to interactive viewing
tm_shape(hdb_risk) +
tm_bubbles()